Large variability in dynamical transitions in biological systems with quenched disorder 
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Coherent oscillatory activity can arise spontaneously as a result of increased coupling in a system 
of excitable and passive cells, each being quiescent in isolation. This can potentially explain the ap- 
pearance of spontaneous rhythmic contractions in the pregnant uterus close to term. We investigate 
the transition to periodic activity using a model system comprising a lattice of excitable cells, each 
being coupled to a variable number of passive cells, whose distribution defines a quenched realization 
(replica) of spatial disorder. Close to the transition between quiescent state and sustained oscilla- 
tions in the system, we observe large fluctuations between different replicas induced by variations in 
the local density of passive cells around an excitable cell. We demonstrate that the disorder-induced 
fluctuations can be described in terms of a simple scaling relation which involves the strength of 
coupling between excitable cells, the mean passive cell density, as well as the logarithm of the system 
size. Our results can be interpreted as suggesting that larger organs will have greater variability in 
the onset of persistent activity. 

PACS numbers: 05.65. +b,87.18.Hf,87.19.R- 



I. INTRODUCTION 

Biological phenomena at different length scales often 
exhibit periodic activity whose frequency can span many 
time scales, and where the rhythmic behavior is crucial 
to the functioning of the relevant systems [1 . For in- 
stance, in the physiological context, oscillations are ob- 
served in variations of intracellular molecular concentra- 
tions [H [3] , activation of pancreatic /3-cell islets that con- 
trol the release of insulin [4 , changes in the activity levels 
of different brain areas [5 , circadian rhythms that con- 
trol the daily sleep-wake cycle [6 and patterns of locomo- 
tion [7 . In a system comprising many different oscillating 
elements, it is crucial for their activity to be synchro- 
nized in order for the generation of rhythmic activity at 
the macro- or systems-level. One of the most significant 
biological contexts in which such synchronization is ob- 
served is in the regular electromechanical contraction of 
the heart. A specialized group of pacemaker cells located 
in the sino-atrial node of the heart coordinates this pe- 
riodic activity in a robust manner [8, 9 . Disruption of 
the coherent collective behavior can result in arrhythmia 
that may be fatal if untreated in many cases [10 . How- 
ever, such clock-like centralized coordinating agencies, al- 
though seen in several contexts [11 , cannot be identified 
for many biological processes, which raises the possibility 
that coherence appears through self-organization among 
the many interacting oscillating elements. 

A biological system whose function is critically depen- 
dent on the occurrence of rhythmic activity is the preg- 
nant uterus, where coherent contractions are observed 
immediately preceding delivery. However, the available 
evidence does not suggest that this is centrally coordi- 
nated by a specialized cluster of pacemaker-like cells, 
as in the heart. For most of pregnancy, the uterus 
does not exhibit any sustained spontaneously generated 



excitation and can be well described as an excitable 
medium [6 . However, on approaching term, episodes of 
periodic activity are observed more and more frequently, 
and they increase in magnitude as well as in duration 
until system-wide powerful contractions eventually ex- 
pel the fetus [12]. As the mechanical activity of the 
uterus is governed by electrical excitation of the tissue, 
the rhythmic contractions are related to the synchroniza- 
tion of oscillations in the transmembrane potentials of the 
cells in the myometrium that forms the bulk of uterine 
wall. In pathological cases (occurring in more than 10% 
of all pregnancies), rhythmic contractions may be initi- 
ated much earlier than normal leading to preterm birth 
[13] which account for a third of infant mortality in the 
USA [14]. It is not yet known why in many cases rhyth- 
mic activity is initiated significantly earlier than normal. 
As there is currently no effective treatment for events 
leading to preterm labor, understanding the dynamical 
processes underlying the onset of self-organized coherent 
rhythmic activity in excitable medium may have poten- 
tial clinical benefits. 

The present work is concerned with the emergence of 
synchronized oscillations in a system of coupled excitable 
and passive cells [T5HT9] that is motivated by the cellu- 
lar organization of the mammalian uterus, where none 
of the constitutive cells can spontaneously oscillate in 
isolation [I2j [131 HO]- Muscle tissue in the uterus, like 
that of many other biological organs, is heterogeneous in 
composition with electrically excitable myocytes (smooth 
muscle cells) forming the bulk but also containing small 
fraction of electrically passive cells like fibroblasts and 
telocytes [131 El EI] • The role of inter-cellular commu- 
nication, including that between dissimilar cell types, in 
generating spontaneous periodic activity in tissue has re- 
cently come into focus |T5l |16] . This is of particular rel- 
evance for the uterus as the electrical coupling between 
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cells through specialized proteins known as gap junctions 
is seen to increase remarkably during pregnancy [22---24^ . 

Our model, inspired by the heterogeneous architecture 
described above, comprises a lattice of diffusively coupled 
excitable cells, each of which may be connected to a num- 
ber of passive cells. As observations made on biological 
tissue do not reveal any regular organization in the con- 
nection density between excitable and passive cells [14|, 
we consider the number of passive cells connected with an 
excitable cell to be randomly distributed about a given 
mean density. Numerical exploration of the model dy- 
namics exhibits a variety of spatio-temporal regimes at 
different values of inter-cellular coupling and average den- 
sity of passive cells [25 , reminiscent of what has been ob- 
served in other systems of coupled electrically excitable 
and passive cells [161 ESI EZ] . More generally, our results 
are relevant for a broad class of systems that comprise 
coupled heterogeneous elements, e.g., in the cardiac con- 
text [28, 29 . In qualitative agreement with biological 
observations of uterine tissue, it is seen from the model 
that strong cellular coupling promotes synchronized os- 
cillations [25 . 

While our earlier work provided a global description 
of the overall dynamics of the model system [25 , us- 
ing a more detailed approach to understand the gene- 
sis of the various dynamical regimes we now point out 
the key role played by disorder. In this paper, we show 
that the local fluctuation in passive cell density can result 
in self-organized emergence of "pacemaker-like regions", 
i.e., clusters of coupled cells that spontaneously oscillate, 
eventually activating the entire system. We stress that 
these are very different from the specialized pacemaker 
cells of the heart whose role cannot be adopted by any 
of the other cells in cardiac tissue under normal circum- 
stances, whereas individual members of the self-organized 
oscillatory groups referred to above are not inherently 
different from the other cells in the system. Thus, in 
the heart, pacemaker function is an intrinsic property of 
certain specialized cells, whereas, in the uterus it is an 
outcome of interactions between heterogeneous cell types. 

A given distribution of passive cells attached to each 
myocyte effectively represents a particular realization of a 
quenched disordered system. The important role played 
by disorder in various phase transitions has been exten- 
sively investigated in physical systems [30l[31]. Here, we 
investigate how the disorder described above influences 
the transition to sustained coherent activity in a vital 
biological organ. Specifically, we consider fluctuations in 
the macroscopic behavior of the system, viz., emergence 
of synchronized oscillations in different replicas, where a 
replica refers to a particular realization of the disorder. 
If excitable and passive cells are coupled with strength 
Cr, the combination is capable of oscillating when the 
number density of passive cells / coupled to an excitable 
cell is within a critical range: fl{Cr) < f < fciCr)- A 
key aspect of the model is that this transition between 
quiescent state and oscillatory activity occurs through 
a subcritical bifurcation ^32l [33], so that the oscillations 



always have a finite amplitude. For f ^ fl^ the system 
exhibits very strong fluctuations from replica to replica as 
the coupling between excitable cells D is varied. This is 
because close to it is not just the mean value of / over 
the entire system but also the spatially fluctuating val- 
ues of local passive cell density that are important. The 
definition of the local density involves a coarse-graining 
length, (i, which is expected to be related to the coupling 
between excited cells as d ex y/D on general physical 
grounds. Therefore, whether a system oscillates for given 
values of Cr and D will depend critically on the details 
of the passive cell distribution in a specific realization. 
More specifically, oscillating groups will emerge from re- 
gions where the local passive cell density is greater than 
The probability of these extreme events, i.e., the oc- 
currence of such regions through statistical fluctuations, 
increases with system size. Based on these arguments 
we derive a scaling relation of the system dynamics as 
a function of system size and cellular coupling strengths 
that we have verified through numerical simulations. 

The paper is organized as follows. In Section II we de- 
scribe the model system followed in Section III by a dis- 
cussion of the properties of an excitable cell coupled to a 
variable number of passive cells. Section IV reports our 
numerical observations of fluctuations in the system be- 
havior across different replicas and explains how regions 
with high local passive cell densities can act as organiz- 
ing centers for oscillatory activity (i.e., "pacemakers") in 
the system. Section V contains the key theoretical result 
of our paper where we derive the probability of occur- 
rence of such pacemaker regions as a function of different 
system parameters. We show that it has a logarithmic 
dependence on system size and obtain a scaling relation. 
Finally, in Section VI, we conclude with a summary of 
our results. 



II. THE MODEL 

As mentioned above, the system under investigation 
comprises electrically excitable as well as passive cells. 
The dynamics of excitable cells can be described by equa- 
tions having the general form 

c„^^ = -/io„(v;,<?i) + re, (1) 

where Ve (mV) represents the transmembrane poten- 
tial, Cm(= 1/iF cm~^) is membrane capacitance density, 
lionifJ^-^ cm~^) represents the total density of currents 
transported across the cell membrane by different ion 
channels, gi are variables describing the gating dynam- 
ics of the ion channels and Fg corresponds to external 
stimulation that could be due to coupling with neighbor- 
ing cells. While the explicit form for the ionic current 
density lion depends on the type of cell being consid- 
ered and the level of biological realism desired, in this 
article we focus on phenomena that occur independent 
of the finer-scale details of specific models and therefore 
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use the generic FitzHugh-Nagumo (FHN) model [9^ . The 
functional form for the ionic current in the FHN system 
is lion = -AVe{Ve-a){l-Ve)^g, whcrc A{= 3) is a pa- 
rameter governing the fast activation kinetics, a{= 0.2) is 
the activation threshold and g represents an effective re- 
polarization current. The time-evolution of g is described 

by 

f = e(K.-5), (2) 

with e(= 0.08) corresponding to a relatively slow rate of 
recovery of the system from the excited state. If the sys- 
tem is capable of oscillation, its period is governed by 
the time-scale ~ 1 /e [9 . For the chosen set of parameter 
values, the dynamics of the FHN system involves either 
relaxation towards a fixed point for a subthreshold per- 
turbation or a finite amplitude excursion before decaying 
to fixed point after an interval (action potential) for a 
strong enough or suprathreshold stimulus. Thus, the ac- 
tive cells are in an excitable regime in which they are not 
capable of exhibiting oscillations spontaneously (i.e., in 
the absence of external stimulation). 

In addition to excitable cells, the system also contains 
passive cells described by the dynamics of the transmem- 
brane potential Vp [34 : 

^=KiVj'-V,) + r,, (3) 

so that activity of the cell tends to decay at an ex- 
ponential rate to its resting value Vj^ (=1-5) with K 
(=0.25) specifying the characteristic time scale of relax- 
ation. Similar to Fg, F^ represents the external stimulus 
to the passive cell. When excitable and passive cells are 
coupled to each other via an electrical conductance (rep- 
resenting gap junctions in biological tissue), the external 
stimulus of each cell corresponds to the input received 
from the other cell so that the respective external stimuli 
terms are: 

Fe =npCr{Vp-Ve), (4a) 

rp=Cr{Ve-Vp). (4b) 

The strength of coupling is represented by Cr while rip 
corresponds to the number of passive cells coupled to an 
excitable cell. 

To investigate spatially extended patterns that are seen 
in biological tissue such as that of the uterus, we con- 
sider a 2-dimensional (2D) system of coupled excitable 
and passive cells. The system comprises a square lattice 
of L X L{= N) excitable cells coupled diffusively to their 
nearest neighbors on the lattice. The coupling results 
in an additional term D\/^Ve in the term Fg in Eq. ([l]) 
where D is the effective diffusion constant (the lattice 
spacing being set to 1) and the Laplacian being numeri- 
cally implemented by a five-point stencil. 

Each excitable cell is coupled to a variable number rip 
passive cells, with each passive cell being connected to a 
single excitable cell. The total number of passive cells in 



the system is M = N x f being the passive cell number 
density. Each passive cell is assigned to a randomly cho- 
sen excitable cell drawn uniformly from the population of 
N cells. We use a binomial distribution for n^, which at 
large TV converges to a Poisson distribution with mean /. 
Each realization of M passive cells randomly distributed 
over N excitable cells, which is an instance of quenched 
disorder in the system, is referred to as a replica. 

We have integrated the system numerically using a 
fourth-order Runge-Kutta scheme, using time step dt = 
0.05. Periodic boundary conditions have been used to 
obtain all the results reported here; our earlier results 
have shown that using no-flux boundary conditions do 
not result in qualitatively different phenomena [25]. Af- 
ter starting from random initial conditions, the system is 
first allowed to evolve for 2 x 10^ time steps (correspond- 
ing to transient behaviors) before recording the data. 



III. GLOBAL OSCILLATIONS AND LARGE 
FLUCTUATIONS 

A. Single excitable cell coupled to passive cells: 
The mean-field limit 

Before proceeding to investigate the spatially extended 
system, we consider the dynamical system comprising an 
excitable cell coupled to one or more passive cells, given 
by Eqs. ([l][4| [l3[l6l|25]. This problem, which formally 
corresponds to zero spatial dimension, can also be viewed 
as the mean-field limit of extremely strong diffusive cou- 
pling between excitable cells such that spatial fluctua- 
tions can be ignored and the lattice effectively behaves 
as a single excitable cell coupled to the average number 
of passive cells, /. Even though an excitable cell couples 
to only an integer number rip of passive cells, the mean 
/ will in general not be an integer. Thus, the integer 
value of rip in Eq. Q can be replaced by a real number 
/, which corresponds to an exact representation of the 
system dynamics in the mean- field limit. 

Fig. [l] (a) shows that the system is capable of exhibit- 
ing spontaneous oscillation over a range of values of Cr 
and /, as has been reported earlier [15]. For a given value 
of coupling Cr, the fixed point solution of the system is 
unstable when /c < / < where the lower and upper 
bounds can be obtained analytically by linear stability 
analysis of Eqs. ([l]|4| around the steady state. A simple 
relation between and Cr is obtained in the limit of 
large i.e., when the passive cell relaxes rapidly. Under 
this condition Eq. ([3| is replaced by an algebraic equa- 
tion such that Fe = fK{Vp^ - Ve)/{1 + K/Cr). Thus, 
as / and Cr appear only via Fg, the critical value of / 
depends on Cr as /^''" = ci^u + ^i,u/Cr^ where ci^ui^i.u 
are independent of Cr- Note that, although the above 
argument is strictly valid only in the large K limit, we 
observe from our numerical results that the dependence 
of the critical values of / on Cr follows the above rela- 
tion even for small K [Fig.[l] (a)]. Thus, fl{Cr) decreases 



4 




FIG. 1: (a) Coupling an excitable cell to / passive cells with strength Cr results in spontaneous oscillations in the region 
bounded between the curves representing fc{Cr) and fciCr)- The shaded region corresponds to the situation where the fixed 
point solution is linearly unstable (analytical result) while the region enclosed by the lines indicate where oscillations are 
observed in numerical simulations. Most of the 2D results shown in this paper are for Cr — 1-9 and / = 0.5, indicated by the 
diamond, slightly below the critical line fc- The representation of the (f,Cr) space shown here demonstrates that the critical 
values of / depend linearly on 1/Cr. The amplitude (b) and frequency (c) of the oscillations vary as a function of / (shown 
for Cr — 1.9), with the region enclosed by the dashed-dotted lines in (b) shown in a magnified view in the inset. Between 
/q = 0.529 and — 0.545 the system exhibits bist ability, which indicates that the Hopf bifurcation at the lower critical value 
of / is subcritical. The dotted line is a schematic representation of the unstable solution. 



monotonically with increasing Cr and approaches a finite 
value as ^ oo. 

While the Hopf bifurcation resulting in loss of stability 
of the fixed point solution at / = f^iCr) is supercritical, 
the bifurcation at / = fl{Cr) is subcritical so that oscil- 
lating solutions of finite amplitude exist for / < f[{Cr)- 
For example, for Cr = 1.9 we observe oscillations when 
f > fi = 0.545 [Fig. [l] (b), inset], while they are never 
observed for f < /q c:^ 0.529. Thus, for mean passive cell 
density in the interval /o < / < Z^, the system is multi- 
stable so that both fixed point and oscillatory solutions 
can be observed numerically [see the narrow interval en- 
closed by the dash-dotted line in Fig. [l] (b)]. 

In this paper we focus on the situation when / < 
/^(Cr). As discussed in detail later for this condition 
we observe large fluctuations between replicas when the 
system is made to undergo dynamical transitions by in- 
creasing the coupling between excitable cells, D. In con- 
trast, when Cr is increased so that / > /^(Cr), increas- 
ing the coupling D results in a much simpler transition 
that approaches the mean-fleld behavior. The replica 
variations are most pronounced when / just exceeds 
limc^^oo /c(^r), the minimal value of /^(Cr). Most of 
the numerical results reported here are obtained with 
/ = 0.5 [Fig.[l](a)], while the minimal value is / = 0.484. 



B. Two-dimensional media 



Biological organs exhibiting electrical activity that are 
often functionally important, such as the heart or the 
uterus, are spatially extended objects, being made of tis- 
sue comprising a large number of excitable and passive 
cells. As mentioned earlier, we have modeled these by 
a 2D system of coupled excitable and passive cells with 
periodic boundaries (Fig. [2|. Fig. [s] shows snapshots of 
activity in two different replicas of such a system as the 
diffusion constant D is increased. We have investigated 
the system for a value of the average passive cell density 
/ = 0.5 and coupling strength between excitable and pas- 
sive cells Cr = 1.9. For these values, f < fl ^ 0.545, so 
that in the mean-field limit the medium does not show 
any oscillation; however, for finite values of I), it is possi- 
ble to observe oscillatory behavior either in local clusters 
or globally as traveling waves (Fig. [3]). We define a non- 
dimensionalized distance fi = {fl — f)/ fl between the 
system under study and the critical system where the 
fixed point becomes unstable in the mean-field limit (for 
/ = 0.5, /i ~ 5.7 X 10~^). For small /i, the behavior 
of the system is sensitively dependent on the particu- 
lar replica chosen (Fig. |3|. This is manifested as large 
fiuct nations in the system behavior, measured by vari- 
ous order parameters defined later. For example, while 
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FIG. 2: Structure of the 2D model system. Excitable cells 
(large circles) lie on a regular square lattice, and are connected 
to their neighbors with a coupling strength D. In addition, 
passive cells (small circles) are connected to excitable cells 
with coupling strength Cr at randomly selected sites of the 
lattice. Each choice of the distribution of passive cells defines 
a replica. 



global synchronization where all cells oscillate with the 
same frequency is seen in one replica [Fig.|3](a), D = 2], 
another replica for the same set of parameters exhibits 
only localized clusters of cells oscillating with different 
frequencies [Fig. [3](b), D = 2]. 

This sensitive dependence of the system dynamics on 
the distribution of Up can also be seen from the {D^ Cr)- 
parameter space diagrams of the 2D system for different 
replicas (Fig. [4|. While for low values of C^, increasing 
D results in complete absence of oscillations in the sys- 
tem ("No Oscillation" or NO phase), for larger values of 
Cr we observe clusters of oscillating cells ( "Cluster Syn- 
chronization" or CS phase) at low with each cluster 
having a characteristic frequency that may differ from 
other clusters. As D is increased, the clusters gradually 
synchronize with each other although isolated regions of 
non-oscillating cells can exist ( "Local synchronization" or 
LS phase), until all cells eventually oscillate at the same 
frequency at a sufficiently high D ("Global synchroniza- 
tion" or GS phase). In numerical simulations, the phases 
were obtained by using suitable order parameters (de- 
fined in the next subsection), and choosing appropriate 
threshold values. As can be seen by comparing Fig.|4](a) 
and (b), the diagrams differ quite significantly in terms 
of the actual dynamical regimes that are observed for the 
same set of values of D and Cr-, illustrating the signifi- 



cant variability from one replica to the other. By aver- 
aging over many such replicas, we can obtain a "mean" 
phase diagram. The diagram obtained here. Fig. [4](c), 
for / = 0.5 is qualitatively similar to the one obtained 
in Ref. [25 for a higher value of / (= 0.7). The region 
corresponding to Coherence or "COH" phase also exists 
in the present case, appearing at large value of D which 
is outside the region of interest of the present paper. 



C. Order parameters 

For a detailed quantitative analysis of the spatiotem- 
poral dynamics of the 2D system, we have used two or- 
der parameters: (i) the fraction of oscillating cells in the 
medium fosc and (ii) the width of the frequency distribu- 
tion as measured by the standard deviation cjy^ where v 
denotes the frequencies of the oscillating cells. The am- 
plitude of an oscillating cell is obtained from the square 
root of the integral of power spectral density (PSD) of the 
corresponding Ve time- series, v being the frequency (> 0) 
at which the PSD is maximum. The fraction of oscillat- 
ing cells in the system fosc = Nqsc/N is the ratio of the 
number of oscillating cells A^osc^ i-e., those having ampli- 
tude higher than a chosen threshold, to the total number 
of cells N . Note that the oscillating cells have approx- 
imately the same amplitude which is a consequence of 
the subcritical nature of the transition across fl (Figjl]). 
This enables a clear distinction between oscillating and 
non-oscillating cells, so that the value of fosc does not 
depend sensitively on the choice of the threshold. 

The different phases shown in Fig. [4] are defined in 
terms of the two order parameters as follows. The LS 
and GS phases both have ^ 0; however, for the for- 
mer fosc < 1 while for the latter fosc — 1- The CS phase 
has a finite cTj^ while the NO phase is characterized by 
fosc — 0. The order parameters for a given system de- 
pend on the exact form of the quenched spatial disorder 
and we can investigate the statistical properties of distri- 
bution of fosc over an ensemble of replicas (Fig. [5]a). At 
very low values of the distribution of fosc is peaked 
about a value close to fosc = /• Upon increasing the 
distribution broadens with a bias towards large values of 
fosc A peak at fosc = 1 develops around D ^ 1. In- 
creasing D further, one observes a higher probability for 
very small values of fosc- At I) ^ 2, the distribution 
has an almost bimodal form, corresponding to realiza- 
tions with either very few oscillating cells, fosc <C 1, or 
almost all cells oscillating fosc ~ 1. At yet larger values 
of the peak at fosc ^ 1 disappears, and the distri- 
bution concentrates around fosc = 0, implying a strong 
reduction in spontaneous activity of the system. This is 
reflected in Fig.[5](b), which shows the ensemble average 
ifosc) (where () denotes an averaging over replicas) and 
the fluctuations about the mean (inset) as a function of 
the coupling D for 2D systems with L = 50 and 100 at 
Cr = 1.9. We observe that at large values of the frac- 
tion of oscillating cells in the system decreases on average 
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(a) D=0.6 D=2.0 D=3.0 V (b) D=0.6 D=1.0 D=2.0 




FIG. 3: Replica-dependent fluctuations close to the transition to sustained oscillations in 2D systems of coupled excitable and 
passive cells, with identical sets of parameters (/ = 0.5, Cr = 1.9, L = 100), and increasing values of the coupling strength D. 
Snapshots of activity Ve for two different replicas are shown in the top rows of (a) and (b). The corresponding time-averaged 
spatial correlation, which provides a characteristic size of the active regions, is shown in the second row. The frequencies of 
individual oscillators in the medium are shown in the pseudocolor plots in the third row (white corresponds to absence of 
oscillation). The last row shows the local density of passive cells averaged over a length scale d indicated below each frame. 
For the replica shown in (a), we observe global synchronization at D = 2 followed by progressive cessation of spontaneous 
periodic activity in the system indicated as a shrinking region of oscillating cells when coupling is increased further. However, 
for the replica in (b), coherent oscillation is not observed as D is increased, with the existing localized oscillating clusters having 
distinct frequencies gradually decreasing in size. 



with and hence is consistent with the mean- field re- 
sult. The standard deviation of fosc becomes extremely 
high in the range 1 < D < 2 which is a signature of the 
large fluctuations between different replicas. Comparing 
between the two curves for L = 50 and 100 shows that 
the larger system has a higher probability of showing os- 
cillations than the smaller one for the same parameter 
values. This dependence on the size of the medium is in 
fact systematic and is discussed in section [V| 



IV. FLUCTUATIONS OF PASSIVE CELL 
DENSITY AND LOCAL ACTIVITY 

The behavior described in the previous section can be 
qualitatively understood by noticing that diffusion effec- 
tively couples excitable cells that are within a distance 
oc \fD of each other. As a consequence, one can expect 



that the behavior of a given excitable cell depends on the 
density of passive cells in its neighborhood characterized 
by the distance y/D. With this motivation, we begin by 
describing the procedure for computing the local density 
of passive cells surrounding a given excitable cell. 



A. Averaging procedure 

We define local passive cell density fid coarse-grained 
over a length scale d as the convolution fid = Kd * rip of 
the spatial distribution of passive cells np{i^j) {i^j = 
1,...,L) over the lattice with a 2D averaging kernel 
Kd{i^ j). For simplicity we have considered separable ker- 
nels, i.e., Kd{i,j) = kd{i/d) x kd{j/d). As the coupling 
term is effectively described by a Laplacian, a natural 
choice for the coarse-graining kernel Kd is the Gaussian 
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FIG. 4: (a-b) Phase diagrams for the 2D system of coupled 
excitable and passive cells with / = 0.5 in two replicas of lin- 
ear dimensions L = 50. Initial values are randomly chosen for 
each point in the diagram. A "mean" phase diagram obtained 
by averaging over many (~ 100) replicas is shown in (c) for 
a system of linear dimension L = 64. The averaging implies 
that if a replica is chosen at random for a given D and Cr, the 
behavior seen in (c) will be observed with high probability. 



kernel Kf given by: 



(5) 



which has been used for the numerical investigation of the 
model system presented in this section. For analytical 
convenience we have also used the simpler "square top 
hat" kernel, K^^^ in section |v| which is based on the 
classical ID top hat filter: 



kJ"{i) = -H{d/2-\A).. 



(6) 



where H{.) is the Heaviside step function, i.e. H{x) = 
when X < and H{x) = 1 otherwise. For the square 
top hat kernel, fidii^j) is the number of passive cells, 
averaged over a square subregion of size Nd = dP centered 
at the point (i, j). In the case of the Gaussian kernel ^ 
the contribution from a site (^^/) to Up at site 
depends on the distance between the two. The variance 
of fid can be written as = a'^/Nd^ where is the 
variance of n^, and Nd (= 4:7r(P) is the effective number 
of sites contributing to fid. 



of passive cells. In addition, the features observed on in- 
creasing the diffusion coefficient D resemble the patterns 
of the local passive cell density seen upon increasing the 
averaging size d. As / < increasing D is expected 
to ultimately suppress oscillation in Fig. [3] in agreement 
with the mean-field analysis. Indeed, we observe that 
at larger values of the number of oscillating cells is 
reduced on average. The coarse-graining procedure re- 
flects this phenomenon as with increased kernel width d 
the variations in fid are reduced significantly. This de- 
creases the probability that a cell will have sufficiently 
high local passive cell density to generate oscillations. 

We define a "pacemaker-like region" to be a group of 
adjacent cells with a local passive cell density fid larger 
than the threshold /^^{Cr)- The oscillatory activity aris- 
ing in these regions may propagate in the form of waves 
to the rest of the system, so that they effectively act as 
pacemakers that are seen in systems having specialized 
coordination centers. 

The number of cells comprising the pacemaker-like re- 
gion varies from replica to replica. When / < we 
expect the size of the region to shrink with increasing 
coarse-graining length eventually disappearing at large 
d (as predicted by mean-field analysis). Thus, for a given 
replica, we can measure the largest value of the coarse- 
graining length, d*, for which a pacemaker-like region still 
exists. This is the smallest value of d for which fid < fc 
everywhere in the 2D system. 

As mentioned above, the coarse-graining length d is 
related to the diffusion coefficient expected to be of 
the form d ~ V^- Therefore, in analogy with d*, we 
can define a value I)* of the diffusion coefficient above 
which the number of oscillating cells in the system goes 
to zero. Fig. [6] demonstrates that the relation between 

and I)* (shown for different replicas and system sizes) 
follows d*^ ^ TD*, where T is the slope of the linear 
fit and defines a characteristic time of the order of the 
slow time scale in the system, ~ 1/e. We note that as 
system size increases, both D* and increase. Thus, for 
systems just below the critical threshold for oscillatory 
activity, i.e., f < f^^, oscillations are seen for a wider 
range of D in larger systems. This size dependence is 
investigated quantitatively in section [V| 



V. STATISTICAL DESCRIPTION 

In this section, we investigate the effect of coarse- 
graining length d on the presence and spatial extent of 
pacemaker-like regions in the 2D system. For analytical 
convenience, we use the 'square top hat' kernel, Eq. (|6|, 
for coarse-graining. 



B. Emergence of pacemaker- like regions through 
diffusion 



A. Distribution of passive cells 



As can be seen from Fig. [Sj large-amplitude oscilla- 
tions occur in regions that have the highest local density 



Consider a 2D lattice with N = excitable cells in 
which a total number of M = f N passive cells are ran- 
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FIG. 5: (a) Probability distribution function of the fraction of oscillating cells constructed from 176 different replicas of a 2D 
system of linear dimension L — 100, with identical values oi f — 0.5 (/i = 5.7 x 10~^), and Cr — 1-9, shown as a function 
of diffusive coupling strength D. (b) Variation of the mean value of /osc, and that of the standard deviation (inset) with D, 
both for systems of size L — 100 (squares) and L = 50 (circles). Starting from the uncoupled case D = 0, as the coupling 
strength becomes larger we observe an increase in the fraction of oscillating cells in the system, in certain cases resulting in 
global synchronization. Further increase of the coupling eventually leads to complete cessation of any activity in the system. 
The fluctuation about the mean is maximum at D ^ 2 (inset), when the distribution exhibits a strong bimodal nature. 



domly distributed. The probability that there are Np 
passive cells in a region containing Nd = (N (0 < < 1) 
excitable cells is given by the binomial distribution: 



M 



c^^(i-c) 



M-N„ 



(7) 



When N ^ oo, this reduces to a Poisson distribution 
with mean fN^ : 



NJ 



For a square top hat kernel, the averaging occurs over 
Nd = d? sites where d is the coarse-graining length and 
the local passive cell density fid = Np/Nd- For a Gaus- 
sian kernel, the corresponding coarse-grained region com- 
prises Nd = 4:71 sites. 



B. Probability of the occurrence of a spontaneously 
oscillating cell in a neighborhood of Nd cells 

Given the probability distribution of passive cells, we 
now determine the probability that a given cell is ca- 
pable of spontaneous oscillations when it is effectively 
coupled through diffusion to a neighborhood consisting 
of Nd cells. This is obtained by considering the probabil- 
ity that the local passive cell density fid = Np/Nd > /c, 



which is given by the cumulative distribution correspond- 
ing to Eq. ([7|: 



M 



Pl= E PNA^p)=kifcNd.M-f^Nd^l)A^) 



where Ix{a^b) is the regularized (incomplete) beta func- 
tion [35] and depends on system size N. In the limit 

function 



and Pi 

, -PjVd' tends to a regularized incomplete gamma 



(9) 



The probabilities obtained for different values of Nd from 
the above expression agrees well with the corresponding 
values numerically obtained from different replicas of the 
2D system in Fig. [7| 

In the limit of large system size N^ can be rewrit- 
ten by expressing the incomplete beta function in Eq. ([8| 
as a continuous fraction [36 and keeping only lower or- 
der terms in N^'^ and N~'^. In addition, neglecting 
( = Nd/N we obtain for any fixed value of /i: 



(10) 



where X = —ja — log(l — /i). This expression can also 
be obtained from Eq. (|9| by expanding the regularized 
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FIG. 6: Relation between the largest values of coarse- 
graining length d* and diffusion coefficient D* for which a 
given 2D system possesses regions with oscillatory activity. 
Each point corresponds to a different replica whereas differ- 
ent symbols represent different system sizes N — . Error 
bars centered around averages of replicas for the same system 
size express the standard deviation. The linear fit between 
d*^ and D*, shown as a broken line, has a slope T 36.6 
that lies in the middle of the range of observed oscillation 
periods (lying between 31 and 40) for different values of the 
diffusive coupling strength in a 2D system with L = 50. 



FIG. 7: Probability P^^ of having local passive cell density 
greater than in a neighborhood of Nd excitable cells as a 
function of Nd. Open (filled) symbols represent numerical ob- 
servations using Gaussian (square top hat) kernel. The broken 
curves indicate the corresponding probabilities obtained from 
the analytical expression of P^^, Eq. Is}. The solid curve 
shows , the probability values obtained for infinitely large 
system, Eq. 



C. Probability of the occurrence of a 
pacemaker- like region in the 2D system 



incomplete gamma function and using the Stirling ap- 
proximation. For small /i, we have A :^ > and on 
introducing a reduced variable X = /j^'^f^^N^ one obtains: 



poo 



-X/2 



(11) 



Note that in order to ensure that all replicas for a given 
value of / have exactly the same mean number of pas- 
sive cells f = Up connected to an excitable cell, we have 
used a binomial distribution of rip for generating the dif- 
ferent replicas in our numerical study. In comparison, 
the Poisson distribution used in our analysis leads to a 
mean fraction of passive cells attached to an excitable 
cell that varies from one replica to another for a given 
value of /. In the limit of large system size, the fluctu- 
ations in the number of passive cells become very small, 
so that both binomial and Poisson distributions lead to 
the same behavior. However, when N is not very large, 
these fluctuations are appreciable when using the Pois- 
son distribution which introduces biases in the results 
(see Figs, [t] and |9|. In addition, for finite A/", corrections 
of order N^/N are expected in the expression involving 
the reduced variable A', Eq. (11). These effects are re- 



sponsible for deviations between the predictions of our 
analysis and the numerical results. 



As mentioned earlier. Fig. [6] shows d*, the largest value 
of the coarse-graining length d for which a pacemaker- like 
region exists in the 2D system, as a function of the diffu- 
sion coefficient. We define the probability II]^^ of having 
a pacemaker-like region in the system as the probability 
of having at least one cell with fid > in a replica of 
size N = cells when the coarse-graining is done over 
a neighborhood of cells. Fig. [S] shows the sigmoid 
form of this probability as a function of Nd^ computed 
by using 400 replicas. To obtain an effective value of d* 
(or N^) for an ensemble of many replicas, we define the 
quantity J* (or TV*) by the condition that =1/2. 

As mentioned earlier, larger systems have higher proba- 
bility of having a pacemaker-region for a fixed value of d 
(or Nd) which is explicitly shown in Fig. [sj 

If we assume that there are r] uncorrelated, i.e. effec- 
tively independent, subsets of size Nd in the 2D system, 
we have: 



l-{l-Pj 



Nd) • 



(12) 



One possible estimate of rj is to assume that the inde- 
pendent subsets are obtained by tiling the system us- 
ing non-overlapping blocks of size Nd^ i.e., rj = N/Nd. 
However, a single pacemaker region may be split among 
several neighboring tiles, which effectively leads to an 
underestimation of the probability H]^^. Another possi- 
ble estimate of r] is to assume that all blocks of size Nd 



10 




FIG. 8: Probability II]^^ of having a pacemaker- like region 
in a 2D system of size A/" = as a function of the coarse- 
graining size Nd. Different symbols represent different sys- 
tems sizes from L = 50 to 6400. The continuous lines indicate 



fits to the theoretical expression for given by Eq. (12) 



FIG. 9: The largest coarse- graining size for which oscillations 
are observed in the system, N^, defined as H^* = 1/2, as a 

d 

function of the system size N. The circles represent numerical 
data shown in Fig. [S] The dashed straight line is a guide 
to the eye indicating the slope l/(A/c). The soli d lin e and 
dashed-dotted line represent the solutions of Eq. ( 13 ) using 
the expressions given by Eq. ^ and Eq. ^ respectively. 



are independent, so that rj = N. However, the high de- 
gree of overlap between neighbouring blocks introduces 
significant correlations between them, leading to an over- 



estimation of n]^^ in Eq. (12). 



In general, we expect that the value of r] will be re- 
lated to N and hy r] = m • N/N^. The existence of a 
number m independent of the system size can be qualita- 
tively understood from the following argument. Consider 
the covariance of the two variables, fd{hj) = ^p{h j)/^d 
and fdij^' ij')! i-e., of the average number of passive cells 
at the two sites (^, j) and {i'^j') separated by a distance 
''')^. It has the exact expression 



I = y'{i-i')^ + {j-f 
ap{l)/Nd where a = (1 — ()f and p{l) 



exp[-{l/2dy] is 



the correlation between the two sites. Thus, how fast the 
two values of fd decorr elate is simply given by To es- 
timate the number of independent subsets containing Nd 
points, we introduce defined as the maximum possible 
correlation between two independent subsets. The corre- 
lation length lo is then defined by p{lo) = so that two 
subsets separated by a distance / > /q are independent. 
This implies the existence of N/Iq independent blocks, 
suggesting in turn that m = Nd/l^. Intuitively, one ex- 
pects lo to be of the order of the width of the Gaussian 
kernel, (i, yielding m ~ 47r. We have numerically ob- 
tained the effective number m by least-square fitting of 
the numerical data shown in Fig. [S] with the theoretical 
expression of Eq. ( 12 ). This gives values of m lying in the 



range 9.5 < m < 16.5, which is in fact consistent with 
the heuristic estimate m = 47r. 



D. System size dependence of the probability of 
occurrence of pacemaker-like region 



We now obtain the cutoff value for the coarse- 
graining size below which oscillations are present in the 
system by solving 



- 1 



N 

(l-P„",)"1 = l/2. 



(13) 



To solve this equation, we assume that the system is large 
enough, so that Pff* can be replaced by P^*, Eq. (10). 



This implies that the only system-size dependence of 

is from the exponent mN/Nd- We thus obtain a linear 

relation between log(A^) and (Fig. [9|: 

li^flm 



/iVX*=21og7V + 21og 



log(2)V2^ 



(14) 



which effectively defines For large systems, the be- 
havior of N2 is well described by the \og{N) term in 
Eq. (14) shown by a straight line in Fig. [9j It is worth 
noting that the log(A^) scaling appears very naturally in 
the general context of extreme value statistics [37] • The 
constant term in the R.H.S. of Eq. (14) is responsible for 
the deviations from \og{N) scaling seen for N ^ 10^. 



VI. SCALING IN THE TRANSITION TO 
ACTIVITY IN SPATIALLY EXTENDED 
SYSTEMS 

The analysis presented in the previous section allows 
us to identify a characteristic coarse-graining size, sug- 
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FIG. 10: (a) Average value and (b) standard deviation of the fraction of oscillating cells in a 2-D system, shown as a function 
of the reduced variable AnDT /N^^ where is defined in Eq. (14). The symbols indicated in the legend correspond to different 
values of L and /, i.e., N and /i. The different curves collapse to the same form (within statistical accuracy) when Nd is large, 
the regime for which the analysis described in the text is valid. 



gesting that the number of pacemaker-like regions in a 
given system is determined by the ratio NjjN^^ where 



is given by Eq. (14). In addition, Fig. 6 shows that 



diffusive coupling D results in coarse-graining of the pas- 
sive cell density over a region of size ^ AttD T^ wit h 
T being the typical oscillation period (subsection |IV B ) . 
Based on this, we expect that the fraction of oscillating 
cells will depend on the reduced variable AttDT/N^. 

Fig. [To] shows the average (a) and the standard devia- 
tion (bjof the fraction of oscillating cells, /osc, averaged 
over many (~ 100) replicas as a function of AttDT/N^^ 
for several values of / (or equivalently, of fi) and L. The 
curves for the mean value of fosc seem to collapse to a 
common form at large where our analysis applies, 
Eq. (10), is not valid when is small; see in particular 



the discussion at the end of subsection VB), supporting 
the conclusion that the properties of the system indeed 
vary with the scaling parameter D/N^. The deviations 
seen in the different curves can be explained as an out- 
come of the large standard deviation in fosc as seen from 
Fig. [To] (b). Remarkably, apart from the curves corre- 
sponding to the average values, we also observe relatively 
good collapse in the curves corresponding to the standard 
deviation of fosc^ at least when D/N^ is not too small. 
This suggests that our analysis may apply not only to 
the averaged value, but also to the second moment of the 
distribution of the fraction of oscillating cells, fosc 

As mentioned in the previous section, in the limit of 
very large system size, is a function of the reduced 
variable ja'^ D / \og{N) . This scaling implies that for a 
fixed value of the coupling between excitable cells, as 



system size is increased and/or the mean value / of the 
passive cell distributions comes closer to the critical value 
there is higher probability of observing sustained os- 
cillatory activity in the system. Alternatively, if / is kept 
fixed at a value smaller than then the larger the sys- 
tem size, the larger the range of values of D over which 
oscillations can be observed. 



VII. CONCLUSION 

In this paper we explore the role of heterogeneity in 
the spatial distribution of passive cells, which introduces 
quenched disorder in a system of coupled excitable and 
passive cells, on the transition between quiescent state 
and oscillatory activity. We begin our analysis with the 
mean-field limit of the spatially extended system, which 
is equivalent to a single excitable cell coupled to a num- 
ber of passive cells. We observe here that the transition 
from quiescent to oscillatory activity is subcritical, i.e., 
oscillations appear at onset with a finite amplitude. We 
next investigate the system at finite values of the cou- 
pling between excitable cells and characterize the coarse- 
graining effect of diffusion. We show that the dynamical 
behavior of the system is related to the local passive cell 
density obtained using a suitable coarse-graining length 
scale. We observe large variability in the dynamical be- 
havior of different replicas, i.e., different realizations of 
the spatial distribution of passive cells. The strong effect 
of this quenched disorder on the dynamics is reminiscent 
of glassy systems. To characterize the properties aver- 
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aged over many realizations of the disorder, we assume 
that the effect of diffusion of strength D is to coarse- 
grain the local passive cell distribution over sites. We 
indeed empirically establish that D oc Nd- From this 
perspective, we analyze the transition between quiescent 
state and sustained oscillations. We obtain a scaling re- 
lation that describes the transition as a function of (i) 
system size, (ii) coupling between excitable cells and (iii) 
average passive cell density. One of the important im- 
plications of this relation is that the occurrence of oscil- 
latory behavior depends on the logarithm of the system 
size N so that increasing TV enhances the probability of 
observing oscillations. 

As mentioned in the introduction, the model system 
that we have analyzed has been motivated by the biolog- 
ical phenomenon of onset of coherent oscillations in the 
pregnant uterus close to term. One of the implications of 
our work is that larger organs may show greater variabil- 
ity in their dynamical behavior for a given set of param- 
eters describing the state of the system. In particular, 



they may be more likely to exhibit oscillations even prior 
to the transition point as a result of spatial fluctuations, 
potentially implying that mammals having bigger uteri 
will be at higher risk of having pre-term rhythmic activ- 
ity. As it is not yet well-understood why in some cases 
periodic dynamical behavior is initiated in uterine tissue 
significantly earlier than normal, our study of the role 
of disorder in creating an effective pacemaker- like region 
giving rise to rhythmic activity in such systems may be 
of potential significance for possible clinical applications. 
Our work also connects the dynamical phenomena seen 
in such biological systems with the study of the role of 
disorder in phase transitions occurring in several physical 
systems, including spin glasses. 
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